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ABSTRACT 

We compute upper limits on the nanohertz-frequency isotropic stochastic gravitational wave background 
(GWB) using the 9-year data release from the North American Nanohertz Observatory for Gravitational Waves 
(NANOGrav) collaboration. Well-tested Bayesian techniques are used to set upper limits for a GWB from su- 
permassive black hole binaries under power law, broken power law, and free spectral coefficient GW spectrum 
models. We place a 95% upper limit on the strain amplitude (at a frequency of yr“*) in the power law model of 
Agw < 1.5 X 10“*^. We hnd that the upper limit on the individual spectral components is consistent with a white 
noise spectrum at frequencies greater than the inverse observation time. For a broken power law model, we 
place priors on the strain amplitude derived from simulations of Sesana (2013) and McWilliams et al. (2014). 

Using Bayesian model selection we find that the data favor a broken power law to a pure power law with odds 
ratios of 22 and 2.2 to one for the McWilliams and Sesana prior models, respectively. The McWilliams model 
is essentially ruled out by the data, and the Sesana model is in tension with the data under the assumption of 
a pure power law. Using the broken power-law analysis we construct posterior distributions on environmental 
factors that drive the binary to the GW-driven regime including the stellar mass density for stellar-scattering, 
mass accretion rate for circumbinary disk interaction, and orbital eccentricity for eccentric binaries, marking 
the hrst time that the shape of the GWB spectrum has been used to make astrophysical inferences. We then 
place the most stringent limits so far on the energy density of relic GWs, Ugw(/)/i^ < 4.2 x 10 yielding a 
limit on the Hubble parameter during inflation of //* = 1.6 x 10“^ mpi, where nipi is the Planck mass. Our limit 
on the cosmic string GWB, < 2.2 x 10“***, translates to a conservative limit of Gp, < 3.3 x 10“*- a 

factor of 4 better than the joint Planck and high-f Cosmic Microwave Background data from other experiments. 
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1. INTRODUCTION 

The 1967 discovery of pulsars (Hewish et al. 1968) inau¬ 
gurated the age of high-energy astrophysics, and in time it 
led to the indirect but compelling confirmation that compact 
binaries lose energy to gravitational waves (GWs) in accor¬ 
dance with general-relativistic predictions (Taylor & Weis- 
berg 1989). It is therefore quite fitting that pulsars should 
now feature prominently in the quest to directly detect low- 
frequency GWs, and to use them as probes of the unseen, 
dynamical Universe. In this article, we describe how one 
can make astrophysical inferences about the underlying GW 
source population using the shape of the nanoHertz GW spec¬ 
trum. Using the 9-year NANOGrav data release, Arzouma¬ 
nian et al. (2015), we provide the first constraints on envi¬ 
ronmental factors which may be contributing the measured 
shape of the spectrum. We expect these novel techniques to 
become standard analysis tools, and as such, provide detailed 
descriptions of our methods. We also provide a summary of 
our results here, so that one may get an overview of our new 
methods and results before delving into the details. 

1 . 1 . Gravitational Wave Detection with Pulsar Timing Arrays 

Sazhin (1978) and Detweiler (1979) realized that GWs 
could manifest as otherwise unexplained residuals in the times 
of arrival (TOAs) of pulsar signals after subtracting a deter¬ 
ministic timing model. This timing model accounts for the 
intrinsic evolution of pulsar spin, radio frequency-dependent 
delays to to interstellar propagation effects, the astrometric 
time delays and advances due to the relative motion of the 
pulsar system with respect to the Earth (indeed, to the ob¬ 
servatory), as well as orbital-kinematic and light-propagation 
effects for pulsars that orbit a binary companion (see, e.g., 
Lorimer & Kramer 2005). Foster &. Backer (1990) pointed 
out that the timing residuals from an array of pulsars (a pulsar 
timing array, or PTA) could be analyzed coherently to sep¬ 
arate GW-induced residuals, which have distinctive correla¬ 
tions among different pulsars (Hellings & Downs 1983), from 
other systematic effects, such as clock errors or delays due to 
light propagation through the interstellar medium. 

Today, three international consortia [NANOGrav 
(McLaughlin 2013), the EPTA (Kramer & Champion 
2013), and the PPTA (Hobbs 2013)] are more than ten years 
into extensive campaigns to search for GWs by timing dozen 
of individual millisecond pulsars (MSPs), in which the best- 
timed have rms residuals less than 100 ns (corresponding to 
GW strain sensitivities ^ 10“^^). The three PTAs collaborate 
and share data under the aegis of the International Pulsar 
Timing Array (Hobbs et al. 2010). 

In order to robustly detect GWs, one must have a thorough 
understanding of the underlying noise in the pulsar timing 
data (see e.g. Cordes 2013; Stinebring 2013, for a detailed 
review). Template matching errors due to radiometer noise 
are uncorrelated in both time and frequency, but pulse-jitter 
noise (Cordes & Shannon 2010) appears to affect all TOAs 
obtained simultaneously in different frequency channels. Cor¬ 
related timing noise with a red power spectrum occurs to vary¬ 
ing degree in different pulsars. Spin noise (Shannon & Cordes 
2010) is achromatic and is much smaller in MSPs compared to 
objects with stronger magnetic fields and longer spin periods. 


Chromatic red noise due to propagation through intervening 
plasmas (ISM, interplanetary medium, and ionosphere) may 
also be present if dispersive delays are not removed perfectly 
or if scattering and refraction effects contribute significantly. 

1 . 2 . The stochastic GW background from SMBHBs 

PTAs are most sensitive to GWs with frequencies on the 
order of the inverse timespan of timing observations, where 
TOA measurement noise averages out most efficiently. The 
strongest expected sources in this band are supermassive 
black hole binaries (SMBHBs) with masses of 1 O*^- 1 O*°M 0 , 
out to z ~ 1 (Rajagopal & Romani 1995; Jaffe & Backer 2003; 
Wyithe & Loeb 2003). The binaries form after the hierar¬ 
chical mergers (Sesana et al. 2004, 2008) of galaxies hosting 
individual SMBHs (as most galaxies are thought to do, cf. 
Kormendy & Ho 2013). The most massive and nearest bi¬ 
naries may be detected individually by PTAs through their 
continuous GW emission. Moreover, the cosmic population 
of SMBHBs may be observed collectively as a stochastic GW 
background composed of the incoherent superposition of sig¬ 
nals from the binaries. Rosado et al. (2015) discuss which 
detection is likely to come first. 

The main focus of this article is the measurement of 
an isotropic stochastic GW background from SMBHBs. 
Anisotropic-background searches based on the formalism and 
techniques developed by Gair et al. (2014); Taylor & Gair 
(2013); Mingarelli et al. (2013) are currently underway and 
will be the subject of a follow-up paper. We briefly con¬ 
sider stochastic backgrounds from relic (or primordial) GWs 
as well as backgrounds from cosmic (super)strings as comple¬ 
mentary studies. Our main focus is to demonstrate how, even 
in the absence of a positive detection, PTA data can be used 
to constrain and characterize the astrophysical processes and 
SMBHB source populations that give rise to the GW back¬ 
ground. 

The simplest characterization of the stochastic GW back¬ 
ground - a power-law Gaussian process with isotropic inter¬ 
pulsar correlations - applies if: 

1 . all binaries are assumed to have circular orbits (so each 
component signal is instantaneously monochromatic); 

2. all binaries evolve through the PTA band due purely 
to GW emission, as opposed to environmental effects 
such as interactions with nearby gas or with stars in the 
galactic nucleus; 

3. all binaries are distributed isotropically across the sky 
in sufficient numbers to fulfill the central limit theorem 
at all frequencies. 

Under these conditions, the observed timing residuals due to 
the GW background are described fully by the (cross-) power 
spectral density 

’y,^, (1) 

where a and b range over the pulsars in the array, 7 = 13/3 for 
a background composed of SMBHBs, and Tab is the Hellings- 
Downs (1983) isotropic correlation coefficient for pulsars a 
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and b (a function of the separation angle between their lines 
of sight, which is normalized so that Taa = 1; see also Min- 
garelli & Sidery 2014). Power-law GW backgrounds are also 
described (independently of observations) in terms of their 
characteristic strain 

/L(/)=Agw(^^) , (2) 

which is related to Eq. (1) by 5,/(/) = TijhdfY/f^) and 
-i = 2-la(a = -ll^ for SMBHBs). 

Recent predictions for the value of Agw, based on mod¬ 
els of SMBH-galaxy coevolution and on observational con¬ 
straints of galaxy assembly and SMBH mass functions, range 
between ~ 10“'^ and 10“'^ (McWilliams et al. 2014; Sesana 
2013; Ravi et al. 2014— hereafter MOP14, S13, and RWS14). 
The MOP14 model assumes that all SMBHBs are in circu¬ 
lar orbits evolving under GW emission alone, as well as a 
single black hole-host correlation from McConnell & Ma 
(2013), yielding an estimate of the stochastic GW background 
that is roughly four times as large as S13 and RWS14. The 
S13 model uses a wide variety of galaxy merger rates and 
empirical black hole-host relations to yield a collection of 
phenomenological SMBHB merger rates, which are used to 
compute a distribution of possible GW signals. Note that 
for this paper we consider distributions from S13 that only 
use the black hole-host galaxy relations of McConnell & Ma 
(2013) and Graham & Scott (2013) The RWS14 model also 
assumes the black-hole-host host correlation of McConnell & 
Ma (2013) but includes the possibility of the SMBHBs evolv¬ 
ing in stellar environments, and accounts for non-zero binary 
eccentricities. Thus, some of these models predict spectral 
densities that deviate from straight power-law behavior at low 
frequencies; in that case, we refer the fiducial Agw to their 
value at a frequency of yr“'. Finally, recent results, (Kor- 
mendy cfe Ho 2013), indicate that the black hole-host corre¬ 
lation’s normalization is being revised to larger values with 
more observations indicating that an even stronger GWB may 
be expected; however, for this work we use the published re¬ 
sults based on McConnell & Ma (2013) to make the most fair 
comparison among models. 

1.3. New results in this paper 

Over the last few years, the three PTAs have reported ever 
improving upper limits on the GW backgrounds of the form 
(1): Agw < 7 X 10“'^ (NANOGrav, Demorest et al. 2013, here¬ 
after DFG13), 6 X 10“'^ (EPTA, van Haasteren et al. 2011), 
3 X 10“'^ (EPTA, Lentati et al. 2015, hereafter LTM15), 
2.4 X 10“'^ (PPTA, Shannon et al. 2013), all quoted at 95% 
confidence and a reference frequency of yr“', although dif¬ 
ferences in the statistical analyses and in the availability and 
selection of pulsar datasets mean that these numbers are not 
entirely homogeneous. It is clear that there is significant ten¬ 
sion between these observational limits and astrophysical ex¬ 
pectations for Agw. It is important to note that a limit on Agw 
does note translate directly to a limit on the SMBHB popula¬ 
tion because of the finite number of pulsars that contribute to 
the limit and the stochasticity of the GW signal itself 

This statement can be made more precise. In Sec. 4.2 of 
this paper we report a new 95% upper limit Agw < 1.5 x 10“'^, 
obtained from the Bayesian analysis of NANOGrav’s 9-year, 
37-pulsar dataset released in 2015 (Arzoumanian et al. 2015). 
This limit is five times more constraining than the same anal¬ 
ysis applied to NANOGrav’s 5-year dataset (DFG13) (see the 


top of Fig. 4 for a comparison of the two posterior probability 
distributions). Now, following Shannon et al. (2013), we can 
assess the consistency of our result with astrophysical GW- 
background models. We find a 0.8% probability that the ob¬ 
served Agw, as characterized probabilistically by its Bayesian 
posterior, is drawn from the amplitude distribution developed 
in MOP14, and a 20% probability that it is drawn from the 
(very similar) RWS14 and S13 distributions. Correspond¬ 
ingly, the two bottom panels of Fig. 4 show that 9-year ob¬ 
servations update significantly the MOP 14 and RWS14/S13 
amplitude priors, much more so than our 5-year dataset. 

In Sec. 4.1 we report also a frequentist, optimal-statistic 
95% upper limit Agw < 1.3 x 10“'^, a fivefold improvement on 
the analogous result of DFG13; however, the optimal statistic 
is problematic in the presence of marginally-detectable GW 
signals, so we offer it only as a proxy for the improving sen¬ 
sitivity of NANOGrav’s observations. 

Stochastic GW backgrounds in the PTA band may also 
originate from quantum fluctuations amplified during infla¬ 
tion (Grishchuk 2005) and from topological broken-symmetry 
remnants such as (Damour & Vilenkin 2001; Olmez et al. 
2010), for which Eq. (1) applies with 7 = 5 (depending on 
the equation of state w) and 16/3, respectively. 

In Sec. 5.2.1 we obtain 95% upper limits Agw < 8.1 x 
10“'® for relic GWs, corresponding to energy-density lim¬ 
its Hgw(/= 1/yr)/!^ < 4.2 x lO”*", where h parametrizes the 
Hubble constant Ho = hx lOOkm/s/Mpc. This limit is a fac¬ 
tor of 2.9 better than limit reported by the EPTA in Lentati 
et al. (2015). We then obtain limits on the Hubble parameter 
during inflation, 77* = 1.6 x 10“^ mpi, where mpi = 1 /s/G is 
the Planck mass, using the method developed by Zhao (201 1). 

In Sec. 5.2.2 cosmic strings, we find Agw < 6 x 10“*® and 
= l/yr)h^ < 2.23 x 10“'®, corresponding to a conser¬ 
vative limit on the string tension of Gp < 3.3 x 10“^ - a 
factor of 4 better than the joint Planck and high-/ Cosmic 
Microwave Background data from other experiments. If we 
then restrict ourselves to a GWB produced by the production 
of large cosmic string loops, as described by Blanco-Pillado 
et al. (2014), then our string tension limit is much more re¬ 
strictive: Gp, < 1.3 X 10“'®, a factor of 6.6 times more con¬ 
straining than an identical analysis performed using the EPTA 
limit. 

1.4. Astrophysical Inference 

The mismatch between observations and expectations can 
be explained in different ways. First, it is possible that the 
astrophysical models and the three assumptions listed above 
are correct, but the background amplitude realized in nature 
lies in the tails of the predicted distributions. This hypothesis 
obviously wanes as upper limits get more stringent. 

Second, it is possible that some of the inputs of the astro- 
physical models are not estimated correctly; we can then use 
GW observations to constrain these inputs. For example, in 
Sec. 5.1.1 we assume measurements of the galaxy mass func¬ 
tion and merger rate, and we constrain the scaling between the 
galaxy bulge mass and the central SMBH mass, which affects 
the observed Agw most significantly, through the distribution 
of binary chirp masses. We find a slight inconsistency be¬ 
tween the scaling relation reported by Kormendy & Ho (2013) 
and our limit, whereas the McConnell & Ma (2013) estimate 
is consistent within its margin of error. The McConnell & Ma 
(2013) black hole-stellar velocity dispersion relation under¬ 
lies the MOP14 predictions, while S13 and RWS14 take into 
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account a variety of alternative black hole-host estimates. 

Third, the simple GW-background characterization that 
yields Eq. (1) may not be realistic, because SMBHBs may 
form with significant eccentricity and retain it into the PTA 
band, distributing GW emission over a range of frequency 
harmonics. Furthermore, environmental effects (interactions 
with stars on centrophilic orbits in galactic nuclei, or with cir- 
cumbinary gas disks) can accelerate the transit of individual 
binary systems through the PTA band (see Secs. 5.1.3 and 
5.1.2 for details and references). These environmental effects 
deplete the GW background at low frequencies where PTA 
measurements are most sensitive (i.e., frequencies ~ the in¬ 
verse observation timespan), so PTA upper limits may yet be 
compatible with the MOP14/S13/RWS14 predictions at the 
higher frequencies where the 7 = 13/3 power law is realized. 

To investigate this point, in Sec. 4.2.2 we reanalyze the 
NANOGrav data using a phenomenological Sijif) in the form 
of an inflected power law, parametrized by a turnover fre¬ 
quency /bend and by a shape parameter k, as proposed by 
Sampson et al. (2015) [see Eq. (24)]. By combining this en¬ 
hanced GW-background model with MOP 14 and S13/RWS14 
amplitude priors, we conclude that the data prefer an inflected 
spectrum to a moderate degree for MOP14, however this 
preference is very weak for S13/RWS14 models. Quantita¬ 
tively, the Bayes factors between enhanced and pure-power- 
law spectral models for each of the two priors are 22 and 
2.2, respectively; graphically, the shading in Fig. 5 represents 
the frequency-by-frequency posterior probability density for 
the GW spectrum, which appears significantly inflected for 
MOP14, and only slightly so for S13/RWS14. The data are 
not sufficiently informative to constrain the amplitude and 
shape of the spectrum jointly in the absence of a compact 
prior, so we cannot produce a unique metric of consistency, as 
we did above in the case of the simple power-law spectrum. 

Beyond this phenomenological characterization, the joint 
(^gw,/bendi tj) posteriors Can be mapped into constraints for 
the SMBHB eccentricities (which we do in Sec. 5.1.3) and for 
the astrophysical variables that govern environmental interac¬ 
tions (Sec. 5.1.2). For the former, we assume for simplicity 
that all binaries had the same eccentricity eo when the semi¬ 
major axis of their orbits was 0.01 pc and that they evolved 
purely by GW emission since; we follow Huerta et al. (2015) 
to construct eccentric binary populations and GW strain spec¬ 
tra. The resulting posteriors on eo indicate that eo > 0.7 
is preferred for the MOP14 prior, and eo > 0.5 is preferred 
for S13/RWS14 (though still consistent with smaller values). 
These limits suggest that either SMBHBs form with rather 
high eo, or that binary eccentricity is not a good explanation 
for the mismatch between Agw observations and predictions. 

To characterize environmental interactions (see Sec. 5.1.2), 
we compute the evolution of orbital frequency due to stel¬ 
lar scattering events and to circumbinary gas disk interac¬ 
tions as d//df cx and respectively, correspond¬ 

ing to K = 10/3 and 7/3 (since for GW-driven evolution 
d//df oc /"/^); the frequency /bend then marks the transition 
between environmentally- and GW-driven evolution. For the 
case of stellar scattering, /bend depends most significantly on 
the mass density p of galactic-core stars. Astrophysical es¬ 
timates for p are quite uncertain with typical values between 
10-10'^ Mq pc“^ Dotti et al. (2007) assuming that a major¬ 
ity of our GW sources come from merging elliptical galax¬ 
ies. Under several simplifying assumptions (e.g., that all bi¬ 
naries have circular orbits, that all galaxies have comparable 


densities, and that only a single environmental effect is ac¬ 
tive), we And that p > lO"^M 0 pc“^ is strongly preferred for 
the MOP 14 amplitude prior, and the data is unconstraining 
for the S13/RWS14 prior. For the circumbinary disk case, 
/bend depends on the accretion rate on to the primary (most 
massive) BH, Mi. The accretion rate of the primary BH is 
a function of Mi, the mass of the primary BH, e, the radia¬ 
tive efficiency parameter with a canonical value of e = 0.1 and 
Kopp is the disk opacity. Mi <x. Mie~^K~^^. Hence Mi takes 
on a range of values, typically 1 O“^M 0 yr“' - 1 Mq yr“', 
see e.g. Di Matteo et al. 2001; Armitage & Natarajan 2002; 
Goicovic et al. 2015. In our analysis, we And that the ac¬ 
cretion rate Mi > lO“'M 0 yr“' is strongly preferred for the 
MOP 14 amplitude prior, and again, the data is unconstraining 
for the S13/RWS14 prior. 

1.5. Plan of the Paper 

This paper presents the first analysis that characterizes the 
spectral amplitude and shape of the GW background; in some 
cases we And tension between observations and predictions 
for the GW background from SMBHBs. When informed with 
astrophysical amplitude priors, the data favor phenomenolog¬ 
ical models that include an inflection point over pure GW- 
driven power laws. We attempt to interpret the phenomeno¬ 
logical posteriors in terms of binary-eccentricity or environ¬ 
mental effects by carrying out our analyses in sequence, in¬ 
vestigating different effects separately. Admittedly, we rely 
on simple analytical models of eccentricity and environments; 
to some extent, this simplicity becomes necessary if we are to 
draw any astrophysical conclusion from data that are not yet 
very informative. Therefore, it will be extremely interesting 
to improve the sophistication of our analysis, and to apply it to 
longer, richer datasets, such as the upcoming NANOGrav 11- 
year dataset, as well as the multiple-PTA datasets assembled 
by the International Pulsar Timing Array. 

The rest of this paper is organized as follows: in Sec. 2 
we describe the NANOGrav 9-year pulsar-timing and dataset; 
in Sec. 3 we discuss our signal and noise models, as well as 
the statistical framework of our analysis; in Sec. 4 we docu¬ 
ment important implementation details, and report our results 
in detail; in Sec. 5 we discuss the astrophysical interpretation 
of our analysis, and derive limits on astrophysical and cosmo¬ 
logical quantities; in Sec. 6 we offer our conclusions. 

2. OBSERVATIONS 

This paper uses the observations from the NANOGrav 9- 
year data release, recently presented in Arzoumanian et al. 
(2015), which contains observations made over a time span 
from 2004 to 2013. Initially the array consisted of 15 
pulsars, and it grew to 37 pulsars over the course of the 
project. The first five years of data on 17 pulsars constituted 
the NANOGrav 5-year data release, previously reported by 
DFG13. In this release, all data have been reprocessed. We 
give a brief overview of the dataset in this section; for details 
we refer the reader to Arzoumanian et al. (2015). 

2 . 1 . Observatories 

The 305 m William E. Gordon Telescope of the Arecibo 
Observatory (AO) was used to observe pulsars with declina¬ 
tions in the range 0° < i5 < 39°; pulsars outside of this range 
were observed with the 100 m Robert C. Byrd Green Bank 
Telescope (GBT) of the National Radio Astronomy Observa¬ 
tory (NRAO). The pulsars PSRs J1713H-0747 and B1937H-21 
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were observed at both AO and the GBT. The typical observa¬ 
tion cadence was about once every month. 

At AO, four receivers were used; 327 MHZ, 430 MHz, 
1400 MHz, and 2100 MHz. Of those, typically two (or more) 
were used in immediate succession within ~1 hour per obser¬ 
vation session, which is possible at AO since the receiver does 
not need to be physically replaced for a receiver change. At 
the GBT, observations were made with two receivers at 800 
MHz and 1400 MHz. Observations were only included in 
the dataset if observations could be made within a time span 
of 14 days; otherwise the observations were discarded for the 
lack of information about variations in the interstellar medium 
dispersion. 

Two sets of nearly-identical data acquisition systems were 
developed specifically for NANOGrav. Early observations 
through 2012.3 (AO) and 2011.0 (GBT) were recorded by 
the Astronomical Signal Processor (ASP) and the Green Bank 
Astronomical Signal Processor (GASP) respectively (Demor- 
est 2007). The later observations beginning at 2012.2 (AO) 
and 2010.2 (GBT) were recorded using the Puerto Rican Ul¬ 
timate Pulsar Processing Instrument (PUPPI) and the Green 
Bank Ultimate Pulsar Processing Instrument (GUPPI) respec¬ 
tively (DuPlain et al. 2008; Ford et al. 2010). These instru¬ 
ments performed real-time coherent dedispersion on the dig¬ 
itized incoming baseband data, and folded the data using a 
pre-computed ephemeris. After RFI excision, polarization 
calibration, and flux calibration, the end product of each in¬ 
strument consisted of total-intensity pulse profiles for a series 
of frequency channels. These profiles were integrated over the 
course of an observation, resulting in one or more subintervals 
of typically 20-30 minutes each. 

2 .2. Time of Arrival data 

TOAs were created using various tools in the PSRCHIVE 
package: the Fourier-domain algorithm of Taylor (1992) to 
calculate TOAs, and denoising the pulse profiles via wavelet 
decomposition. Offsets resulting from latencies between dif¬ 
ferent observing systems were determined as overall timing- 
model-fit-parameters, which were taken into account when 
doing the timing noise analysis. 

For each pulsar, TOAs were calculated for all frequency 
channels recorded from a given receiver. The effect of 
time-varying dispersion was taken into account by including 
“DMX” parameters in the timing model (DFG13), which es¬ 
sentially allows for an extra delay proportional to 1 jv^ to be 
fit for, where v is radio frequency. Pulse shape evolution with 
frequency is taken into account differently than in DFG13. 
Instead of including a phase offset per frequency channel, a 
heuristic mitigation procedure is used that parameterizes the 
profile evolution with frequency, the details of which are in 
Arzoumanian et al. (2015). 

In summary, similar to the NANOGrav 5-year release, the 
NANOGrav 9-year release consists of high-quality, publicly 
available' TOAs for 37 pulsars, produced per frequency sub¬ 
channel. 


which follows that of Arzoumanian et al. (2015). We con¬ 
tinue our discussion with the Bayesian approach in Sec. 3.2, 
and the frequentist approach in the form of the optimal cross- 
correlation statistic in Sec. 3.3. 


3.1. Likelihood 

We start our discussion by decomposing our Ntoa timing 
residual data for a single pulsar (5t in its individual constituents 
as follows: 

(5t = Me-i-Fa-i-t/j-i-n. (3) 

The term Me describes inaccuracies in the subtraction of the 
timing model, where M is called the timing model design ma¬ 
trix, and e is a vector of small timing model parameter off¬ 
sets. The term Fa describes all low-frequency signals, in¬ 
cluding low-frequency (“red”) noise, with a limited number 
of Fourier coefficients a. Our harmonics are chosen as inte¬ 
ger multiples of the harmonic base frequency 1 /T, with T the 
length of our dataset (either of a single pulsar, or the entire 
array of pulsars). The matrix F then has alternating sine and 
cosine functions. We note that this is just a particular choice 
of rank-reduced basis, and we could have chosen many others 
without influencing our results. The term t/j describes noise 
that is fully correlated (correlation coefficient of 1) for simul¬ 
taneous observations at different observing frequencies, but 
fully uncorrelated in time. The matrix U is an A^toa x Aepoch 
matrix that maps the Aepoch observation sessions to the Atoa 
TOA s. The vector j describes the white noise per observation 
session that is fully correlated across all observing frequen¬ 
cies. The last term, n, describes Gaussian white noise that is 
assumed to be left in the data, after correcting for all known 
systematics. 

The white noise is assumed to be an uncorrelated, Gaussian 
noise process, with variance-covariance: 

(n,n,) = ^ Nij^k = Y. + QlSij), (4) 

k k 


where Ek and Qk are the TEMPO and TEMP02 “EFAC” 
and “EQUAD” parameters for each observing backend k, Sjj 
is the Kronecker delta function, and W = diagfa?} are the 
TOA uncertainties. The matrix elements (i,j) apply to the 
TOAs corresponding to the observing system labeled by k. In 
practice we cannot fully separate various contributions to our 
TOAs, so we have to take into account that our corrections 
for the various terms of Eq. (3) are inaccurate. We do this by 
constructing our likelihood from the noise-mitigated timing 
residuals: 

r = St-Me-Fa-Uj, (5) 


where r is our best approximation of n, given our knowledge 
of all the noise and signal parameters. The likelihood can now 
be written down as: 


p{6t\e,a,i,(j)) 


exp(-ir^A-'r) 
^/det{2^TN) ’ 


(6) 


3. DATA ANALYSIS METHODS 

All the data analysis methods we use in this manuscript 
(both the Bayesian and frequentist methods) are effectively 
canied out in the time-domain. We start in Sec. 3.1 by 
defining the likelihood function and introducing our notation. 


http :/ /data . nanograv . org 


where N = J2k^k represents the total effect of all white noise. 
We have collectively denoted all parameters not directly rep¬ 
resented by e, a, and j as f. Henceforth, we shall refer to e 
as the hyperparameters. We group the reduced-rank signals as 
follows: 


T=[M F U], 


b = 


e 

a 

J. 


1 


(7) 
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which allows us to elegantly place a Gaussian prior on the 
coefficients of these random processes. The prior covariance 
is: 

"oo 0 O' 


B = 


0 if 0 , 

0 0 J_ 


( 8 ) 


resulting in a prior: 


p{b\cf)) = 


exp(-jb^B 'b) 

y/det(2TTB) ’ 


(9) 


where oo is a diagonal matrix of infinities, which effectively 
means we have a uniform unconstrained prior on the tim¬ 
ing model parameters e. As described in Arzoumanian et al. 
(2015), this representation allows us to analytically marginal¬ 
ize Eq. ( 6 ) times Eq. (9) over the waveform coefficients b 
of the noise, resulting in a drastic reduction of the dimension¬ 
ality of our posterior (Lentati et al. 2013; van Haasteren & 
Vallisneri 2014, 2015): 


piSt\cf)) = 


exp(-jSt^C 'i5t) 

^det(27rC) ’ 


( 10 ) 


with C = N+TBT^. The Woodbury matrix identity (Wood¬ 
bury 1950) can be used to evaluate Eq. (10) efficiently. 

The parameters that describe B are the hyperparameters 4>. 
The hyperparameters of the diagonal matrix S are the per- 
backend TEMP02 ECORR parameters. The matrix ip repre¬ 
sents the spectrum of the low-frequency noise and the stochas¬ 
tic gravitational waves, and it therefore contains terms corre¬ 
lated between pulsars. Denoting pulsar number with (a,b), 
and frequency bin with we can write: 


“ ^abPi^ij Pai^ab^ij ■; ( 11 ) 

where is the spectrum of the low-frequency noise of pulsar 
a, p is the spectrum of the GWB, and Tab is the signal corre¬ 
lation matrix. The elements of the signal correlation matrix 
consist of the overlap reduction function for a GWB signal, 
which is a dimensionless function that quantifies the corre¬ 
lated response of the pulsars to a stochastic GWB (Mingarelli 
& Sidery 2014). The quantities p, and rjai can be either mod¬ 
eled as independent model parameters (i.e., per frequency), 
or as modeled spectral density with a specific shape (e.g., a 
power-law model). We note that p can be represented by 
a block-diagonal matrix, where each block corresponds to a 
specific frequency bin; all frequencies are theoretically inde¬ 
pendent degrees of freedom. 


3.2. Bayesian analysis 

As an alternative to the orthodox frequentist approach to 
data analysis, Bayesian inference is a method of statistical 
inference in which Bayes rule of conditional probabilities is 
used to update one’s knowledge as observations are acquired. 
Given a model H, model parameters 0, and observations T), 
we write Bayes rule as: 

Pr{e\VM)Pr{V\H)=Pr{V\e,H)Pr{Q\H ), (12) 

where Pr (Q\'D,'H) = Pr{Q) is the posterior (probability dis¬ 
tribution) of the parameters, Pr ?() =L(0) is the likeli¬ 

hood, Pr(0|'H) = 7r(0) is the prior (probability distribution), 
and Pr (D|'H) = Z is the marginal likelihood or evidence. 

The left-hand side of Eq. (12) can be regarded as the “out¬ 
put” of the Bayesian analysis, and the right-hand side is the 


“input”. Indeed, provided we have a generative model of our 
observations (meaning we can simulate data, given the model 
parameters), we know the likelihood and prior. However, for 
parameter estimation we would like to know the posterior, and 
for model selection we need the evidence. 

Eor parameter estimation, the evidence Z is usually ig¬ 
nored, and one can use L( 0 ) 7 r( 0 ) directly to estimate con¬ 
fidence intervals. Typically one provides confidence intervals 
for single components and pairs of elements of 0. This in¬ 
volves an integral over Pr{Q) over all but one or two param¬ 
eters, a process called marginalization. When 0 is higher¬ 
dimensional, Monte-Carlo sampling methods are typically 
used to perform this multi-dimensional integral. We use 
Markov Chain Monte Carlo methods in this work to sample 
the posterior distribution. 

Model selection between two models "Ho and Hi can be 
carried by calculating the “Bayes factor”: the ratio between 
the evidence for the two models. Assuming we have a prior 
degree of belief of how likely the two model are (Pr(Ho) and 
Pr{Hi)), we can write: 


Pr(Hi\V) 

Pr{Ho\V) 


= BioiV) 


PHHi) 

PriHoY 


(13) 


where BioCD) = ZijZo is the Bayes factor, and O is the odds 
ratio. The odds ratio can be obtained by calculating the evi¬ 
dence Z for each model separately (e.g. with Nested Sam¬ 
pling or thermodynamic integration), or by calculating the 
Bayes factor B between two models directly (e.g. with trans- 
dimensional markov chain Monte Carlo methods). 


3.3. Optimal cross-correlation statistic 

The optimal statistic (Anholm et al. 2009; Demorest et al. 
2013; Chamberlin et al. 2015) is a frequentist estimator of the 
isotropic GWB strain-spectrum amplitude in the weak-signal 
regime. The estimator maximizes the likelihood of Eq. (10), 
and it can be written as: 


2 _ 


— 

gw 






(14) 


where = (i5ta(5tj) is auto-covariance matrix of the residu¬ 
als of a single pulsar. This is C of Eq. (10) for a noise-only 
model. The term A^^Sab = Sah = {StaStJ) represents the sig¬ 
nal covariance between pulsar a and pulsar b. As described 
in the previous section, our model contains no other signals 
with a non-zero correlation coefficient between different pul¬ 
sars. The normalization of the optimal statistic is such that 

Eollowing Chamberlin et al. (2015), we also quote expres¬ 
sions for the variance, and the signal-to-noise ratio (S/N, in 
power) of the optimal statistic. TTie variance of the optimal 
statistic in the absence of a GWB signal is given by: 



-1 




(15) 


Although this expression is not valid in general, in the weak- 
signal regime, in which the cross correlated power is ignored 
from the variance, this can be used as an approximation of 
the variance in The S/N for a given signal and noise 















NANOGrav Nine-year Isotropic GWB Limit 


7 


realization is; 


CTo 




(j^ab 


tr 




which has an expectation value of 


(p) =^gw E 


1/2 ’ 


1/2 


tr 




ab 


(16) 


(17) 


These expressions generalize the detection significance esti¬ 
mator provided in Jenet et al. (2005), properly taking into 
account the spectrum of the signal and the noise, as well 
as details such as the irregular sampling. The S/N here has 
the same meaning, which if interpreted as a zero-mean unit- 
variance normal distribution, can be used to place upper-limits 
on the GWB amplitude. Clearly, this interpretation is only ap¬ 
propriate in the weak-signal regime, but it serves as an inde¬ 
pendent sanity check for our other methods. Setting our statis¬ 
tical significance threshold a = 0.95, we can place a one-sided 
upper-limit as: 


A2i=i2^ + V2aoeifc-'[2(l-a)]. (18) 

We note that all the usual caveats of frequentist-type upper- 
limits apply to this methodology as well, as no prior informa¬ 
tion is used. For instance, it is possible to set an upper-limit 
of Aui = 0 in a dataset without a (detectable) signal, which 
theoretically happens with probability a. 

It was shown in Chamberlin et al. (2015) that the optimal 
statistic is identical to the cross-correlation statistic used by 
DFG13. This alternative interpretation of the optimal statistic 
allows us to obtain a measure of the cross power between pul¬ 
sars. The cross power is the amount of correlated power be¬ 
tween the timing residuals of different pulsars, and one would 
expect this cross power to follow the Hellings and Downs 
cross-correlation signature for a detectable GWB signal. The 
cross power and the uncertainty estimates are given by 


^ab ~ 


St^PAPb'^t, 


tr 




'^0,ab - (tr 




- 1/2 


(19) 

which is independent of any specific overlap reduction func¬ 
tion and A^^FabSab = Sab- One can then fit these correlation 
coefficients assuming a particular overlap reduction function 
by minimizing 


x^=E 


ab 


Cab 


ab 


<^0,ab 


( 20 ) 


where Tab are the cross-correlation coefficients given by the 
overlap reduction function for an isotropic GWB. It can easily 
be shown that the best fit value of Agw is then the optimal 
statistic value of Eq. (14). 


4. RESULTS 

4.1. Optimal Statistic Analysis 

The optimal cross correlation statistic was applied to the 
full 9-year NANOGrav dataset. The maximum-likelihood 
single-pulsar noise values were obtained by independent noise 
analyses on each pulsar. The maximum-likelihood amplitude 



C [deg] 


Figure 1. (Top Panel): Cross-correlated power vs. angular separation ob¬ 
tained through an Optimal Statistic Analysis. The dashed red curve shows 
the maximum-likelihood amplitude mapped onto the Hellings and Downs 
coefficients. (Bottom Panel): Histogram of the number of pulsars in each 
bin (red right axis) and a weighted histogram (blue left axis) using the l-cr 
uncertainties in the cross correlation as weights. 


and SNR of this search are Ag* = 8.9 x 10“'® and p = 0.8, 
respectively, indicating little evidence for the expected GWB 
cross correlations. The resulting upper limit using this method 
is Ag®''® = 1.3 X 10“*® at a reference frequency of yr“*, which 
is 5.4 and 1.5 times more stringent than the limits using this 
method presented in DFG13 and LTM15, respectively. The 
corresponding signal-to-noise ratio is 0.8, indicating little ev¬ 
idence for the expected cross correlations from a GWB. It 
should be noted that the limiting technique presented in An- 
holm et al. (2009) and DFG13 does not strictly have proper 
frequentist coverage in the presence of any measurable GWB 
signal^; therefore, this limit will serve as a proxy to our im¬ 
proved sensitivity and not a strict upper limit. 

As a by-product of the optimal statistic analysis, we can 
obtain the maximum-likelihood cross-correlation values for 
each pulsar pair in the analysis. In Figure 1 we plot the cross- 
correlated power vs. angular separation in the top panel. We 
have binned the values into 10 degree bins using a weighted 
averaging technique. The red curve shows the maximum- 
likelihood correlated power fit. It is clear that the cross- 
correlated power is still consistent with zero signal and that 
we have much more sensitivity at values of C < 100°, which 
is a by-product of the fact that our best timed pulsars all lie at 
a similar position on the sky. This is illustrated in the bottom 
panel of Figure 1, where we plot the histogram of the number 
of pulsars in each ten-degree bin (red) as well as the most sig¬ 
nificant bins (blue). To calculate the most significant bins we 
use the one-sigma uncertainties on each maximum-likelihood 
cross-correlation value as weights in the histogram. 

^ In the limit of zero GWB signal, this method actually over covers in 
the frequentist sense. We have found hy comparison with other frequentist 
bounding techniques and Bayesian upper limits that the 5-year limit pub¬ 
lished in DFG13 is indeed robust and does not suffer to under coverage. 
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4.2. Bayesian Analysis 

We have used the Bayesian data analysis techniques de¬ 
scribed in Section 3 to place upper limits on the GWB param¬ 
eterized using the standard power law, spectral components 
and a broken-power-law. Bayesian data analysis of modern 
PTA data sets is quite difficult to perform in general due to 
the large parameter spaces (~ 500 for the current NANOGrav 
dataset assuming a power law red noise and GWB spectrum) 
and likelihood evaluation time. Even with the efficient meth¬ 
ods described in Section 3 we are still prohibited from carry¬ 
ing out a full Bayesian search where we vary all noise param¬ 
eters and GW parameters for all pulsars simultaneously. 

In this work we ameliorate these two problems by using 
a two-tiered approach to noise modeling and GWB analysis. 
To avoid including all noise parameters in our GWB analysis, 
we first carry out single-pulsar noise analyses on each pul¬ 
sar in which we include EFAC, EQUAD, and ECORR for 
each backend/receiver system and red noise parameterized 
as a power law. We then perform the GWB analysis while 
holding all EEAC, EQUAD, and ECORR values fixed to their 
mean value from the single-pulsar analysis but allowing the 
red noise and GWB parameters to vary. By holding these 
white noise parameters fixed, we reduce the computational 
burden since large matrix products can be pre-computed. Eur- 
thermore, by holding all white noise values fixed we reduce 
the number of free parameters drastically from ~ 500 to 
2A^psr + -^gwj where Apsr and Agw are the number of pulsars 
in the array and number of parameters describing the GWB, 
respectively. 

We further reduce the computational burden in two ways: 
by ignoring cross correlations and by only using a subset of 
pulsars. We choose to ignore cross correlations in this work 
as they do not contribute to the upper limit in the absence of 
a signal and only serve to greatly increase our computational 
burden by requiring that we invert a large dense matrix for ev¬ 
ery iteration in the analysis. That the cross correlations have 
no bearing on the upper limit has been argued in Shannon 
et al. (2013) and further shown in LTM15. A real GWB sig¬ 
nal will reveal itself as a strong common red noise term well 
before the cross correlations become detectable, and since we 
see no evidence for a common red noise term, we are justified 
in dropping these terms. Lastly, we choose to only include a 
subset of pulsars in our analysis as not all pulsars contribute 
to our upper limit either due to short timing baselines or large 
measurement errors. To choose this subset of pulsars, we first 
carry out our single pulsar analyses mentioned above but now 
include an extra red noise process with power spectral index 
fixed to 13/3 (i.e., SMBHBs). The pulsars are then sorted 
based on their single pulsar upper limits on the GWB. We 
then carry out the GWB analysis mentioned above to com¬ 
pute an upper limit using an increasing number of pulsars in 
the sorted list. This process is continued until the upper limit 
saturates. In other words, including more pulsars beyond this 
point does not change the upper limit. Using this method we 
choose to use the 18 pulsars shown in Table 1 . 

To compute the posterior probability of Eq. (10) and to 
map out the multi-dimensional parameter space we make use 
of the pulsar timing data analysis suite PAL2^ in conjunc¬ 
tion with the parallel-tempering Markov Chain Monte-Carlo 
(PTMCMC) code"^. The details of the PTMCMC sampler can 

^ https : //github . com/jellisl8/PAL2 

^ https : //github . com/jellislS/PTMCMCSampler 


Table 1 

Pulsars used in GWB analysis. 


PSR Name 

[10-‘5] 

RMS“ [fis] 

Timing Baseline [yr] 

J1713+0747 

1.96 

0.116 

8.76 

J1909-3744 

4.5 

0.081 

9.04 

J1640+2224 

11.8 

0.158 

8.9 

J1600-3053 

12.3 

0.197 

5.97 

J2317+1439 

13.6 

0.267 

8.87 

J1918-0642 

16.0 

0.34 

9.01 

J1744-1134 

16.1 

0.334 

9.21 

J0030+0451 

19.8 

0.723 

8.76 

J0613-0200 

23.4 

0.592 

8.58 

J1614-2230 

23.9 

0.189 

5.09 

B1855+09 

26.6 

1.338 

8.86 

J1853+1303 

31.0 

0.235 

5.6 

J2145-0750 

33.0 

0.37 

9.07 

J1455-3330 

37.9 

0.694 

9.21 

J1012+5307 

43.3 

1.197 

9.21 

J1741+1351 

56.8 

0.103 

4.24 

J2010-1323 

83.3 

0.312 

4.08 

J1024-0719 

92.9 

0.28 

4.01 


“Weighted root-mean-square of epoch-averaged post-fit timing residuals. 
See Arzoumanian et al. (2015) for more details. 

be found in Appendix C of Arzoumanian et al. (2014). The 
parameters and prior ranges used in the analysis are shown 
in Table 2. As shown in the table, we have chosen to use 
uniform priors on both the red noise and GWB amplitude pa¬ 
rameters. The GWB amplitude prior is chosen so that it is 
proper at Agw = 0, that is to say, it must converge to a finite 
value in the limit of zero amplitude, otherwise the upper limit 
would depend on the lower bound of the prior which is un¬ 
desirable. However, for the red noise amplitude we have no 
such restriction. A log-uniform prior would result in the most 
unbiased parameter estimation and a uniform prior would bias 
the parameter estimation towards higher red noise amplitude 
and lower red noise spectral indices (since there is a strong 
covariance between these two parameters); however, the log 
uniform prior would cause intrinsic red noise to be modeled 
by the GWB amplitude which is also undesirable. Because 
of these considerations we choose uniform priors on both as 
it gives equal weight to both red noise and the GWB. We 
also find that this prior results in a much more robust upper 
limit when using different numbers of frequencies in the rank- 
reduced approximations of Eq. (9). 

4.2.1. Power-law and Spectral Limits 

When computing upper limits on the dimensionless strain 
amplitude we fix the spectral index (13/3 in the case of 
SMBMBs) and adopt a uniform prior on Agw. When perform¬ 
ing a spectral analysis we again use priors that correspond to 
a uniform prior on Agw, this results in a prior that is uniform 
in the square root of the power spectrum coefficients. 

Figure 2 show the results of the power law and spectral 
analysis along with relevant astrophysical model predictions. 
The solid black and long dashed black lines are the 95% up¬ 
per limits from the spectral and power-law analyses, respec¬ 
tively. The green, blue, and red shaded regions are the one- 
sigma prediction on the strain spectra from MOP 14, RWS14, 
and S13, respectively. We find an upper limit on the dimen¬ 
sionless strain amplitude of Agw* = 1.5 x 10“'^, a factor of 
1.6 better than the most constraining published upper limit to 
date (Shannon et al. 2013), a factor of 2 more constraining 
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Table 2 

Summary model parameters and prior ranges. 


Parameter 

Description 

Prior 

Comments 

White Noise 

Ek 

EFAC per backend/receiver system 

uniform in 10,10] 

Only used in single pulsar analysis 

Qk 

EQUAD per backend/receiver system 

uniform in logarithm 1-8.5,—5] 

Only used in single pulsar analysis 

Jk 

ECORR per backend/receiver system 

uniform in logarithm 1-8.5,—5] 

Only used in single pulsar analysis 

Red Noise 

'^red 

Red noise power law amplitude 

uniform in [10-20,10-“] 

1 parameter per pulsar 

'Yred 

Red noise power law spectral index 

uniform in 10,7] 

1 parameter per pulsar 

GWB 

Agw 

GWB power law amplitude 

uniform in [10-'^10^“] 

1 parameter for PTA for power-law models 

7gw 

GWB power law spectral index 

delta function 

Fixed to different values depending on analysis 

pi 

GWB power spectrum coefficients at frequency i/T 

uniform in [10“'®, 10“®]“ 

1 parameter per frequency 

A 

GWB broken power-law amplitude 

log-normal^ for models A(B) 


K 

GWB broken power-law low-frequency spectral index 

M-14.4(-15), 0.26(0.22)) 
uniform in 10,7] 

1 parameter for PTA for broken power law models 

1 parameter for PTA for broken power law models 

/bend 

GWB broken power-law bend frequency 

uniform in logarithm 1-9,—7]^^ 

1 parameter for PTA for broken power law models 


“ The prior uniform in is chosen to be consistent with a uniform prior in Ag„ for the power law model since ip,- oc Ag, 


* These values are quoted in log base 10 and are obtained from MOP14 and S13. 

We choose different prior values on /bend when mapping to astrophysical model parameters as described in Section 5. 



Figure 2. Strain amplitude vs. GW frequency. The solid black and long 
dashed black lines are the 95% upper limits from our spectral and power-law 
analyses. The red, blue and green shaded regions are the one-sigma predic¬ 
tions from the models of S13, RWS14, and MOP14. The green shaded region 
uses the simulation results from MOP14, but replaces the fit to the GWB pre¬ 
dictions used in that paper with the functional form given by Eq. (24). 

than the recent EPTA upper limit of LTM15, and a factor of 
5 more constraining than the 5 year data release upper limit 
when applying the same Bayesian analysis. Furthermore, we 
find a slightly less constraining upper limit when using the 
free spectrum model (power-law equivalent upper limit of 
2 X 10“'^). This is to be expected since the free spectrum 
model has many more degrees of freedom (we use 50 free 
amplitudes for each of the 50 frequencies in this case) over 
the power law parameterization (1 degree of freedom). Thus, 
since the power law model can leverage extra information at 
all frequencies, as opposed to the spectrum model where each 
frequency is independent of the others, more constraining up¬ 
per limits are expected from a power law model. We also find 
that the upper limit on the strain spectrum from the spectrum 
analysis is consistent with white noise (i.e., /!*’'‘‘®(/) cx /^^^) at 


frequencies > 3/T, where T is the length of the longest set of 
residuals in the data set, which indicates that our GWB upper 
limits are coming from the three lowest frequency bins. This 
behavior is to be expected since we have several well timed 
pulsars that do not span the full 9-year baseline (see Table 1) 
and thus will have peak sensitivity at frequencies greater than 

i/r. 

From inspection of Figure 2 we see that our 95% upper limit 
is within at least the 2-sigma confidence region of all three as¬ 
trophysical models and is sensitive to a potential turnover in 
the spectrum due to environmental coupling factors. We wish 
to determine the level of consistency between our data and 
the power-law models displayed in Figure 2. To accomplish 
this we follow the method applied in Shannon et al. (2013). 
Given that we have a model M for the value of the GW ampli¬ 
tude Agw whose probability distribution function is denoted 
/>(Agw|M) and that we have a probability distribution func¬ 
tion for Agw given the data, denoted p{Ag^\d), where d repre¬ 
sents the data, the probability that we measure a value of Agw 
greater than that predicted by the model, A^, is given by the 
law of total probablility 

/ OO pOQ 

/7(Agw|M)^fAgw / />(Agw|^f)iiAgw. 

OO 

( 21 ) 

Therefore, low values of P{Ag^ > A^) indicate that the range 
of Agw that is consistent with our data is inconsistent with 
the model M, and vice versa. To carry out this procedure 
the distribution p{Ag^\d) is simply the marginalized poste¬ 
rior distribution when using the uniform prior on Agw. We 
use log-normal distributions to model the MOP14, S13, and 
RWS14, models. Since the models of RWS14, and S13 pre¬ 
dict nearly the same GWB amplitude distribution (assuming 
a power-law only) we make no distinction between these two 
models. Furthermore, the model distributions on Agw, given 
by log-normal distributions have mean and standard devia¬ 
tions of (-14.4,-15) and (0.26,0.22) for the MOP14 (here¬ 
after model A) and S13/RWS14 (hereafter model B) mod- 
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Figure 3. Upper limit on the GWB as a function of power spectral index. 

els, respectively. Using the aforementioned distributions and 
Eq. (21) we find that our data are 0.8% and 20% consistent 
with models A and B, respectively, under the assumption of a 
power-law. This indicates that either the assumptions that go 
into these models are incorrect, our universe is a realization 
of the GWB that has an amplitude in the tail of the probability 
distributions mentioned above, or that environmental effects 
are depleting SMBHB sources at low frequencies making the 
power-law assumption faulty. The implications of this last 
point are discussed in Section 5. 

In addition to our power law limits on the stochastic back¬ 
ground (i.e., strain spectral index -2/3), we have also com¬ 
puted the upper limit on a the GWB for a range of differ¬ 
ent spectral indices. In Figure 3 we plot the upper limits 
obtained at varying spectral indices (red points) vs. power 
spectral index. We also provide the best fit model for the up¬ 
per limit as a function of power spectral index where we find 
cx cx This differs fromDFG13 where they 

find (X r“, arguing that this is due to the fact that the sen¬ 
sitivity is dominated by the lowest frequency of l/T. Our fit, 
giving a slightly weaker dependence on a is consistent with 
what we have seen above, namely that our limits are not com¬ 
pletely dominated by the lowest frequency. 

In a Bayesian analysis, the posterior distribution is the prior 
distribution updated by the data. Here we illustrate this by 
comparing our power-law upper limits, using identical meth¬ 
ods, on the 5-year (DFG13) and 9-year (Arzoumanian et al. 
2015) NANOGrav data releases. The results of this compar¬ 
ison are shown in Figure 4 where we plot the marginalized 
posterior distributions of logjoAg* for the 9- and 5-year data 
releases in blue and red, respectively. The gaussian prior dis¬ 
tributions described above are shown in green for model A 
and model B. For the uniform prior case we see quite a dra¬ 
matic improvement (i.e., the factor of 5 mentioned above) in 
the upper limits. For model A; the 5-year dataset does some¬ 
what inform the prior, whereas the 9-year data set results in 
a posterior that is largely inconsistent with the prior distribu¬ 
tion. For model B the 5-year data set do not inform the prior at 
all, whereas the 9-year data set does indeed update the prior. 

4.2.2. Broken Power-law Limits 

We place constraints on the strength of environmental cou¬ 
pling effects that will likely affect our GWB signal at low 
frequencies (i.e., large orbital separations) via a simple pa¬ 
rameterization of the GWB spectrum that allows for a “bend” 
frequency at which there is a transition from environmentally- 
driven evolution to GW-driven evolution. The following dis- 
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Figure 4. Marginalized posterior probability density of logjoAgw computed 
using the nine (blue) and five (red) year NANOGrav data releases for uniform, 
MOP14 model gaussian, and S13/RWS14 model gaussian prior distributions. 
The gaussian priors are shown in green. 


cussion and analysis techniques are based on Sampson et al. 
(2015). Here we give a brief overview of this more general¬ 
ized GWB spectrum. 

The characteristic amplitude of a stochastic background 
from an ensemble of SMBHBs in circular orbits is (Phinney 
2001; Sesana et al. 2008; McWilliams et al. 2014) 


hciff = 



(PN dt 
dzdM dt d\n f 


hHf), 


( 22 ) 


where d^N/(dzdAidt) is the differential number of inspi- 
raling binaries per unit z, M and f, where z is the redshift, 
Ai = (miOT 2 )^/^/(mi -I-OT 2 )^/^ is the chirp mass of the binary, 
and t is the time measured in the binary rest frame. The 
dt/d\n f term describes the frequency evolution of the binary 
system, and h(f) is the strain spectrum emitted by a single 
circular binary with orbital frequency //2. Typically, it has 
been assumed that the binary is purely GW-driven which re¬ 
sults in our usual expression for hdf) given in Eq. (2); how¬ 
ever, physical mechanisms other than GW radiation that are 
necessary to drive the binary to coalescence (Milosavljevic 
& Merritt 2003) will change the frequency dependence (i.e., 
the dt/d\nf term) of this equation (see Colpi 2014, for a 
review of SMBHB coalescence). Following Sampson et al. 
(2015) we can generalize the frequency dependence of the 
strain spectrum to 


dt 

d\nf 


= f 




(23) 


where i ranges over many physical processes that are driv¬ 
ing the binary to coalescence. If we restrict this sum to GW- 
driven evolution and an unspecified physical process then the 
strain spectrum is now 


hc{f)=A 


(f/fyrT 

(l+(/be„d//)")‘^"’ 


(24) 
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where /bend and k are the parameters that encode information 
about the physical processes (other than GW radiation) driv¬ 
ing the binary evolution. As mentioned above, there could be 
many physical processes contributing to the frequency evolu¬ 
tion of the SMBHB system; however, at cmi'ent sensitivity it 
is very unlikely that our data can distinguish them. Thus we 
follow Sampson et al. (2015) and adopt this slightly simpli¬ 
fied spectrum to place constraints on possible environmental 
coupling mechanisms. 

The above discussion has focused on the frequency evolu¬ 
tion of SMBHB. The other piece of the equation is the merger 
rate of SMBHs, which will set the overall amplitude scale for 
the GWB somewhat independently (assuming the last parsec 
problem is solved) of the physical mechanisms that drive the 
system to merger. Here we use the same log-normal distri¬ 
butions introduced in Section 4.2.1 for models A and B as 
prior probability distributions for the GWB amplitude A in 
Eq. (24). In the following analysis we also fix a = -2/3 and 
use uniform priors on logjo/bend G [-9,-7] and k S [0,7] un¬ 
less stated otherwise. 

Here, we have run an identical analysis to that of Section 
4.2.1 except that the GWB spectrum model is now that of 
Equation (24) and we adopt the aforementioned priors on A, 
/bend, and K. Eigures 5 and 6 show the results of this anal¬ 
ysis. Eigure 5 shows the posterior probability density of the 
GWB spectrum defined in Equation (24) with model Aon the 
left and model B on the right. This probability density is con¬ 
structed by drawing values of A, /bend, and k from the joint 
probability distribution shown in Eigure 6, constructing the 
spectrum at each frequency via Equation (24) and then pro¬ 
ducing a histogram of the spectral power at each frequency. 
The solid black lines in Eigure 5 represent the 95% credible 
region and the median of the GWB spectrum. The dashed line 
is the upper limit on Agw using the purely GW-driven spec¬ 
trum (i.e., no transition frequency) and the gaussian ampli¬ 
tude priors from models A and B, respectively. Lastly the thin 
solid black line is the 95% upper limit on the GWB spectrum 
from the spectral analysis presented in Section 4.2.1 . By in¬ 
specting the inferred GWB spectrum one can determine that 
the data prefer a GWB spectrum that has a definitive transition 
from GW-driven to environmentally-driven within the pulsar 
timing frequency band for Model A, whereas the data does not 
significantly constrain the shape of the spectrum for model B. 

This can be seen further in the joint posterior distributions 
in Eigure 6 in which the probability distributions (blue) for the 
bend frequency parameter, /bend, and spectral index parame¬ 
ter, K, are significantly different from the prior distribution 
(green) for model A and not significantly different for model 
B. When this same analysis is carried out on the 5-year data 
release we find that the data can only slightly inform the prior 
on /bend for model A and gives no information on the other pa¬ 
rameters in either model A or B. This, again, indicates that the 
9-year data release provides us with much more information 
about the shape of the GWB strain spectrum. 

Einally, we can be more quantitative and apply Bayesian 
model selection to this problem by computing the Bayes fac¬ 
tor between the broken-power-law (Equation (24)) and pure- 
power-law (Equation (2)) parameterizations for both models 
A and B. When this analysis is carried out we anive at Bayes 
factors of 22.2 ± 1.1 and 2.23 ±0.15 for models A and B, 
respectively. These Bayes factors were computed using par¬ 
allel tempering and a custom thermodynamic integration im¬ 
plementation (See Sec. 6.1.2 of Cornish & Littenberg 2015). 

This analysis shows, for the first time, that PTAs are en¬ 


tering a regime where even in the case of a non-detection 
meaningful constraints can be placed on the dynamical history 
of the SMBHB population. Eurthermore, this analysis shows 
that when placing upper limits to make statements about the 
full range of astrophysical merger scenarios it is no longer 
valid to consider only the classic strain amplitude, but one 
must instead frame the question in terms of measuring the 
amplitude and shape of the GWB spectrum. As we have seen 
in the above analysis and as can be seen clearly in Eigure 6, 
given a model for the SMBHB merger physics (i.e., a prior on 
A) and discarding the assumption of a purely GW-driven sig¬ 
nal (i.e., a broken-power-law model), it is very difficult to rule 
out any of the GWB amplitude parameter space with any cer¬ 
tainty unless one has strong a-priori knowledge on the shape 
of the spectrum. However, we can begin to place constraints 
on the environmental coupling effects that drive the system to 
the GW-dominated regime. 

5. DISCUSSION 

5.1. Astrophysical Model Limits 

While the parameter estimation of the previous section is 
interesting in and of itself, some of the most interesting sci¬ 
ence available from the NANOGrav data is accessible only 
by relating these parameters to properties of the source pop¬ 
ulations. Here we attempt to interpret the phenomenological 
posteriors on the GW spectrum in terms of black hole-host 
galaxy relations, environmental effects, or binary-eccentricity 
by carrying out our analyses in sequence, investigating dif¬ 
ferent effects separately. In Section 5.1.1 we use the results 
of our power-law analyses of Section 4.2.1 to provide con¬ 
straints on the parameters of scaling relations between host 
galaxies and black holes. We then go beyond the assumption 
of a power-law spectrum in Section 5.1.2 to investigate how 
our constraints on the shape of the characteristic strain spec¬ 
trum from Section 4.2.2 map to constraints on the environ¬ 
ment of SMBH binaries. Einally in Section 5.1.3, we probe 
the eccentricity of binaries before they entered the PTA band. 

5.1.1. Constraints on Host Galaxy - Black Hole Scaling Relations 

If the gravitational wave spectrum is assumed to be created 
by an ensemble of binary SMBHs that are formed following 
galaxy mergers (spectral index of-2/3), and whose evolution 
is assumed to be dominated by GW emission throughout the 
PTA-sensitive waveband, then we can trace the expected bi¬ 
nary SMBH population using observations of galaxy merger 
rates, the galaxy stellar mass function, and the black hole-host 
galaxy relation. This is the approach taken in SI 3, and Ravi 
et al. (2015). Assuming equal fractional uncertainties in these 
parameters, the black hole-host galaxy relation will have the 
largest impact on the predicted level of the GW background. 
This is due to the much stronger dependence of the GW back¬ 
ground on the chirp mass of each source than on the number 
of sources. 

Middleton et al. (2015) shows that it is difficult to extract 
information from PTA limits without making significant as¬ 
sumptions about the shape of the black hole merger rate den¬ 
sity. If instead a galaxy merger rate density calculated from 
observed galaxy parameters is used as a proxy for the black 
hole merger rate density, then a limit on the GW background 
can be directly translated into a limit on the input galaxy pa¬ 
rameter spaces (Simon & Burke-Spolaor 2015). Eor this pa¬ 
per, we focus specifically on the scaling relation between host 
galaxy properties and black hole mass (e. g. Haring & Rix 
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Figure 5. Probability density plots of the recovered GWB spectra for models A and B using the broken-power-law model pai'ameterized by (Agw, fhend, and k) 
as discussed in the text. The thick black lines indicate the 95% credible region and median of the GWB spectrum. The dashed line shows the 95% upper limit 
on the amplitude of purely GW-driven spectrum using the Gaussian priors on the amplitude from models A and B, respectively. The thin black curve shows the 
95% upper limit on the GWB spectrum from the spectral analysis. 



Figure 6. One- and two-dimensional posterior probability density plots of the spectrum model parameters Agw, /bend, and k. In the one-dimensional plots, we 
show the posterior probability from the 9-yeai‘ data set (blue), the 5-year dataset (dashed red) and the prior distribution used in both analyses (green). In the two 
dimensional plots we show a heat map along with the one (solid), two (dashed), and three (dash-dotted) sigma credible regions, model A is on the left and model 
B is on the right. 


(2004), McConnell & Ma (2013)) as it is the observed pa¬ 
rameter that is most easily constrained by NANOGrav data. 
Specihcally, we constrain the M, -Mbuige relation; 

logioM, = a-i-/31ogio(Mi„/gs/lO"M0) . (25) 

In addition to a and jS, observational measurements of this 
relation also fit for e, the intrinsic scatter of individual galaxy 
measurements around the common a, /? trend line. In prac¬ 
tice, a and e have the greatest impact on predictions of Ag*, 
and all observational measurements agree with « 1. 


PTAs are most sensitive to binary SMBHs where both black 
holes are ^IO^Mq (e.g. Sesana et al. (2008)). Therefore 
M ,-Mtuige relations that are derived including the most mas¬ 
sive systems are the most relevant to understanding the pop¬ 
ulation in the PTA band. Several recent measurements of 
the M ,-Mhuige relation specihcally include high-galaxy-mass 
measurements, e.g. those from Brightest Cluster Galaxies 
(BCGs). As these hts include the high-mass black holes that 
we expect to dominate the PTA signals, we take these as the 
“gold standard" for comparison with PTA limits (Kormendy 
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Figure 8. The above plot shows the translation of the 95% upper limit on 
Agw, Fig. 4, into the parameter space a-e, which characterizes the black 
hole-host galaxy relation as described in Equation. (25). The parameter space 
above the line is inconsistent with the power-law analysis of the S13 model, 
as described in Simon & Burke-Spolaor (2015). Observational measurements 
of this parameter space are shown with errorbars. 


Figure 7. The above plot shows the translation of the marginalized posterior 
distribution of Agw, Fig. 4, into the black hole-host galaxy parameter space, 
which is characterized by an intercept a, a slope 0, and an intrinsic scatter e. 
Flat priors are used for a, 0, and e. 0 is not informed by the distribution of 
Agw, while both a and e are, with a limit on a being more strongly set. The 
curves show the 1, 2, and 3 (j contours. Relevant observational measurements 
are also shown, with McConnell & Ma (2013) in blue and Kormendy & Ho 
(2013) in magenta. Since 0 is not strongly informed by the upper limit, we 
can set an upper limit in a-e space by marginalizing over 0. That upper limit 
is shown in Fig. 8. 

& Ho 2013; McConnell & Ma 2013, e.g.). 

The translation of an upper limit on Agw to the black hole- 
host galaxy parameter space is calculated as follows: 

p(a,/3,e|PTA)cx JdO p(9) piA,^(a, (26) 

where the posterior of Agw, p (Agw(Q;,/3,e,0)|PTA), is the 
marginalized posterior distribution of Agw, which is shown in 
Fig. 4; Agw (a,P,e, 9) is the prediction of Agw calculated from 
models similar to SI 3; 9 represents the galaxy stellar mass 
function and the galaxy merger rate; and p (a,/3,e|PTA) is 
the marginalized posterior distribution of the black hole-host 
galaxy relation, which is shown in Fig. 7. For this analysis, 
we use two leading measurements of the galaxy stellar mass 
function, Ilbert et al. (2013) and Tomczak et al. (2014), and 
two measurements of the galaxy merger rate, Robotham et al. 
(2014) and Keenan et al. (2014), as the basis for simulating a 
local population of binary SMBHs. A flat prior is used for a, 
P, and e, and the posterior on Agw using a uniform prior, as 
seen in Fig. 4, is directly translated into this parameter space. 
The result of which is shown in Fig. 7. /3 is clearly not in¬ 
formed by a PTA posterior, but the combination of a and e 
are, with the strongest limit being set on a. 

Fig. 8 shows the translation of our posterior on Agw into a-e 
parameter space with observational measurements of the pa¬ 
rameters from Kormendy & Ho (2013) and McConnell & Ma 
(2013). Assuming a power-law analysis of the S13 model, as 
described in Simon & Burke-Spolaor (2015), there is a slight 
inconsistency between our upper limit and the Kormendy & 
Ho (2013) measurement. Potential solutions to an inconsis¬ 
tency include; a significant number of black hole binaries do 



Figure 9. By introducing the parameter Tstaii, as described in Simon & 
Burke-Spolaor (2015), we can start to explore the inconsistency of our up¬ 
per limit with power-law models for the GW background. In the above plot, 
we allow Tstaii to vary while using the M* -Mbuige relation Kormendy & Ho 
(2013). The probability of is a direct translation of the posterior on Agw 
from Fig. 4. The blue line is the 95% lower limit on Ltall, which we set 
at 0.73 Gyr. While this is not sufficiently constraining to make meaningful 
astrophysical statements, this pai'ameter may be useful for future PTA upper 
limits. 

not reach the GW-dominant regime in our assumed timescale 
(they “stall"); the ‘classical’ assumption of a power-law strain 
spectrum in the PTA band is incorrect and in fact there is a 
turn-over in the strain spectrum at lower frequencies (see Sec. 
4.2.2); or that the measured astronomical parameters are not 
correct for the population of binary SMBHs in the PTA band. 

As the possibility for a different strain spectrum curve 
is discussed in Sec. 4.2.2, let us explore the potential for 
‘stalling’ within the model described so far in this section. 
Using the galaxy merger rate density as a proxy for the black 
hole merger rate density implies as assumption that the events 
occur at a similar cosmological time. If there was signifi¬ 
cant stalling in the binary black hole population, then these 
events would be offset in cosmological time by some ‘stalling 
timescale’. There is then nothing in the model to allow for 
anything other than an efficient solution to the ‘last-parsec’ 
problem Begelman et al. (1980). As done in Simon & Burke- 
Spolaor (2015), we introduce a variable, Tstaii, which is a mea- 
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sure of this offset in time between the assumed galaxy merger 
rate density and the black hole merger rate density used to 
model a GW background. Fig. 9 shows the translation of the 
posterior distribution of Ag„, Fig. 4, into a probability dis¬ 
tribution of Fstaii using the Kormendy & Ho (2013) measure¬ 
ments of the M,-Mbuige relation. Using this we set a 95% 
lower limit on Fstaii of 0.73 Gyr, which is not sufficiently con¬ 
straining to indicate which of the many suggested ‘solutions’ 
to the ‘last-parsec’ problem are not occurring for the systems 
in the PTA band. However, this parameter may be useful for 
future PTA upper limits as a probe of the level of ‘stalling’ in 
the binary black hole population. 

5.1.2. Constraints on binary environmental influences 

The cores of galactic merger remnants can harbor stars 
with little angular momentum and almost radial trajectories 
which intersect the central galactic region (centrophilic or¬ 
bits). These stars can undergo three-body interactions with 
the resident supermassive black-hole binary, causing the stars 
to be ejected, which results in energy and angular momentum 
being extracted from the black hole system, and leading to bi¬ 
nary hardening (Quinlan 1996).^ Additionally, the formation 
of circumbinary gaseous disks can lead to interactions which 
extract energy and angular momentum from the binary orbit, 
driving it towards smaller orbital separations (Ivanov et al. 
1999). We expect that, in the type of gas-poor galaxies which 
dominate the nanohertz GW signal, hardening from stellar 
scattering will dominate over circumbinary interactions, but 
we consider both in the following. We begin with a discus¬ 
sion of how these environmental mechanisms drive the evolu¬ 
tion of the binary, then provide constraints on the frequency 
at which the characteristic strain spectrum exhibits a turnover 
from the analysis in Sec. 4.2.2. We finish by mapping these 
frequencies to constraints on the astrophysical environment of 
SMBH binaries emitting GWs in the nanohertz band. 


as (Sesana 2013) 


da 

dt 


2Mi 

IJ- 


(aao) 


1/2 

5 


(28) 


where Mi is the accretion rate onto the primary black hole, /r 
is the binary reduced mass, and qq is a characteristic orbital 
separation at which the enclosed disk mass equals the mass 
of the secondary black hole. The latter can be expressed as 
(Ivanov et al. 1999) 



where a is a disk viscosity parameter. Mi 2 are the binary 
black hole masses; Me = 4-ttGMi /ckj is the Eddington ac¬ 
cretion rate of the primary (kj is the Thompson opacity coef¬ 
ficient); and rg = 2GMi/c^ is the Schwarzschild radius of the 
primary. 

As in the stellar scattering case, we can rearrange this equa¬ 
tion to determine the orbital frequency evolution. In this 
model, d//df (x and so k = 1 13 for a-disk binary in¬ 
teractions. 


Constraints on spectral turnover frequency — We define the 
spectral turnover frequency in the obvious way to mean the 
frequency at which the characteristic strain spectrum exhibits 
a change in slope. If the low-frequency slope is positive, this 
will correspond to the point at which the spectrum is maxi¬ 
mized. One must note that setting / = /bend in Eq- (24) does 
not maximize the characteristic strain spectrum. Rather the 
turnover frequency will be a function of both /bend and k: 


Environmentally-driven orbital evolution — We use the formal¬ 
ism of Quinlan (1996) to define the evolution of a (circular) 
binary due to three-body stellar scattering events, where 


d /1\ _ GpH 
dt \a) a ’ 


(27) 


where a is the binary orbital semi-major axis, p is the mass 
density of galactic core stars, H is a dimensionless hardening 
rate which takes a value of ^ 15, and tr is the velocity disper¬ 
sion of core stars. Using Kepler’s third law, we can rearrange 
this equation to solve for the rate of orbital frequency evolu¬ 
tion, which gives d//df cx /E^. Since the binary orbital evolu¬ 
tion will be due to a combination of environmental influences 
and GW emission, df/dt is actually a sum over all mecha¬ 
nisms, as in Eq. (23). We know that the rate of frequency 
evolution due to GW emission is d//df oc /'E^ (Peters & 
Mathews 1963), hence, in the language of the parametrized 
spectrum model in this paper given in Eq. (24), k= 10/3 for 
binary hardening by three-body stellar scattering. 

Likewise, the evolution of a circular binary due to circumbi¬ 
nary disk interaction is modeled within the a-disk formal¬ 
ism (Ivanov et al. 1999; Haiman et al. 2009; Sesana 2013) 


^ We assume that all galactic merger remnants maintain the same mass 
density of core stars throughout the binary merger. The subtleties of loss- 
cone replenishing impact the evolution of the binary and of the central density 
profile within a factor of ~ 2 (a few at most), as shown by Sesana & Khan 
(2015) and Vasiliev et al. (2015). 


\ 

/turn=/bend( ^-ij • (30) 

We can combine our measurements of /bend and k from the 
analysis in Sec. 4.2.2 to compute the probability distribution 
of spectral turnover frequencies. We find that placing numer¬ 
ical constraints on /bend is difficult as the posterior is heavily 
influenced by the prior, namely the upper and lower bounds 
for the uniform priors used in this analysis. In the following 
analysis we set the lower bound on /bend by requiring that the 
power at / = 1 /Tjpan differs by no more than 10% from a pure- 
power-law for any value of k. By using this prior, we ensure 
that we can recover a pure-power law spectrum (in our fre¬ 
quency range) for any value of k. The upper bound of /bend is 
chosen based on the specific environmental coupling mecha¬ 
nism we are considering. 

We emphasize that the probability distributions of /bend, k, 
and /turn are not distributions of these parameters over the 
SMBH binary population. Rather our posterior distribution 
is illustrating the spread of our beliefs in the measurement of 
the single /bend and k model parameters. 

Constraints on environmental parameters — We can now extract 
astrophysical constraints from our constraints on the transi¬ 
tion and spectral turnover frequencies. By equating the rate 
of orbital evolution, da/dt, due to environmental mechanisms 
and GW emission, we can deduce the characteristic transition 
frequency, /bend, between these influences. We firstly consider 
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stellar-scattering, for which the transition frequency is given 
by 

/bend = 3.13 nHz X (31) 

= 3.25 nHz x (32) 

where M is the total binary mass. Mg = M/(10* Mq), 
q = M2/M1, qr = q/(l+q)^, CT200 = cr/(200 km/s), pg = 
p/(10^ M0pc“^), and His = H/15. In the second line we 
have used the M-a relation (McConnell et al. 2011) to re¬ 
place velocity dispersion, where Mg « 1.9 <j\qq. It follows 
from the McConnell et al. (2011) M-a relation that the ve¬ 
locity dispersion term on the first line of Eq. (31) has the 

mass scaling a^Q^ oc , which is a small modifica- 

—2/5 

tion to the exponent of the other mass term. Mg ' . Hence, 
/bend a relatively weak dependence on the mass scaling 
of a. Nevertheless we can express the transition frequency 
in terms of variables of a parametrized M-a relation, where 
log[o(M/MQ) = fl-l-l?log[0(7200, such thatMg = lO^'^CTjog, and 

/bend = 3.13 nHz X 

(33) 

Finally, the weak scaling of /bend with H, and the < 20% devi¬ 
ations of this parameter away from 15 seen in numerical scat¬ 
tering experiments (Sesana & Khan 2015) justifies our keep¬ 
ing this parameter fixed at its fiducial value of His = 1 ■ Astro- 
physical estimates on p are quite uncertain with estimated val¬ 
ues around 10-10"^ Mq pc“^ for typical environments (Dotti 
et al. 2007). The variation of estimates over several orders 
of magnitude is why we choose to investigate only p in our 
stellar-scattering constraints. 

The equivalent transition frequency for a-disk interaction 
is 

/bend « 0.144 nHz (34) 

where aq is as previously defined. We adopt the fiducial value 
of the disk viscosity parameter a ~ 10“^ used in Ivanov et al. 
(1999). The very weak dependence of the bend frequency on 
this parameter, /bend oc will significantly dampen the in¬ 
fluence of any deviations from this fixed value. Hence, in our 
constraints on the influence of circumbinary disk interactions, 
we only vary the accretion rate of gas onto the primary black- 
hole, Ml, of which estimates in the literature vary over several 
orders of magnitude - typically 1O“^M0 yr“' - 1 Mq yr“\ 
see e.g. Di Matteo et al. 2001; Armitage & Natarajan 2002; 
Goicovic et al. 2015. 

Equations (31)-(34) indicate the GW frequency at which a 
single circular binary will transition between being environ¬ 
mentally driven and being GW driven. Of course, /bend is not 
the quantity that can be extracted from a spectral analysis - 
/turn is what we can measure. Our analysis of the stochastic 
GW background has provided us with constraints on the char¬ 
acteristic transition frequency of the entire population, and 
thus the turnover of the spectrum. Hence, our path to pro¬ 
viding constraints on environmental parameters requires us to 
numerically construct characteristic strain spectra for popu¬ 
lations of SMBH binaries in contact with their environment. 
We can then construct numerical mappings between the envi¬ 
ronmental parameters of interest (core stellar mass density, p, 
and primary black hole accretion rate. Mi) and the turnover 
of the spectrum. We use the formalism of Phinney (2001) via 


Eq. (22), where the differential comoving number density of 
merging binaries per redshift and component masses is con¬ 
structed as in McWilliams et al. (2014): 

^^^^ocP(z)/(Mi)/(M 2), (35) 

where P(z) encapsulates redshift dependent factors, and 
(j>{M^ 12 }) is the number density of black holes of a certain 
mass, which is given by a (redshift dependent) Schechter 
function modified at the high-mass end by a lognormal com¬ 
ponent to accommodate recent high mass BCG (brightest 
cluster galaxy) discoveries (Lin et al. 2010). 

The combined influence on the binary orbital evolution of 
GW emission and environmental couplings is modeled with 
the sum in Eq. (23), where either stellar scattering or disk 
interactions are included i.e. we consider one mechanism at 
a time. By considering all binary environments to have the 
same p or Mi, we can determine the spectral turnovers from 
the numerically computed strain spectra, iterating over many 
values of these environmental parameters to deduce a map¬ 
ping. We draw binary systems from the ranges z S [0,1], 
Ml G [10^,10^°] Mq, and q G [0.1,1]. The results for our 
fiducial assumptions are shown in the top panel of Fig. 10, 
with the stellar density required to give a certain turnover fre¬ 
quency shown in the left panel, and the primary accretion rate 
required to give a certain turnover frequency shown in the 
right panel. 

In the lower panels of Fig. 10 we plot the posterior dis¬ 
tributions of the stellar density, p, for stellar hardening, and 
mass accretion rate. Mi for circumbinary disk interaction. In 
this analysis we perform the Bayesian parameter estimation 
for fixed values of k corresponding to the appropriate values 
for stellar hardening (k = 10/3) and circumbinary disk inter¬ 
action {k, = 7/3). These posteriors are constructed by first 
converting the marginalized distributions on /bend to a dis¬ 
tribution for /turn via Eq. (30) and then using the empirical 
mapping to convert /um to the appropriate astrophysical pa¬ 
rameter. Again, we do not place numerical confidence limits 
on p or Ml since the data does not constrain the prior distribu¬ 
tion at large values. Nonetheless, from inspection of Fig. 10 
we see that the MOP14 model heavily prefers p // 10“* M0pc“^ 
and Ml > 10“' M0yr“', while the S13 model is largely uncon¬ 
straining for both mechanisms. Typical densities of massive 
elliptical galaxies at the MBH influence radius is ~ 10^ Mq 
pc“^, making the MOP14 model hard to reconcile with ob¬ 
servations, even if we consider that massive ellipticals were 
likely factor 2-3 more compact at z = 1 (Sesana, unpublished). 
Our results approach the upper end (for the MOP 14 prior) 
of the expected range of M, IQT^Mq yr“' - 1 Mq yr“', ob¬ 
served in the local universe and predicted via simulations, 
see e.g. Di Matteo et al. 2001; Armitage & Natarajan 2002; 
Goicovic et al. 2015. Furthermore, Dotti et al. (2015) predict 
that Ml 10“' Mq yr“' for BH masses of 10® Mq and red- 
shifts z < 1; however, these are average accretion rates, and 
short, episodic accretion triggered by galaxy mergers could 
occur at a higher rate. 

We go beyond the fiducial assumptions for the case of stel¬ 
lar hardening since it is the most likely environmental cou¬ 
pling mechanism for SMBH binaries. When we increase the 
low-mass cutoff of systems which contribute to the charac¬ 
teristic strain budget this further constrains the stellar mass 
density. This is seen most easily in Eqn. (31), where one must 
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Figure 10. {top)\ Empirical mapping from /tum to p {left) and My (right), (bottom)'. Posterior distributions for the mass density of stai's in the galactic core 
(left) and the accretion rate of the primary black hole from a circumbinary disk (right). These distributions are constructed by first converting the marginalized 
distribution of /bend to a distribution of /um via Eq. (30), and then using the empirical mapping described in the text to convert from /urn to the astrophysical 
quantities p and M \, respectively. 


raise the stellar mass density to match a corresponding in¬ 
crease in binary mass so that the transition frequency is main¬ 
tained. Furthermore, modeling the distribution of black holes 
masses in Eq. (35) without the lognormal component or red- 
shift evolution will increase the contribution of lower mass 
binaries to the GW strain budget, leading to smaller stellar 
mass density constraints than reported in Fig. 10. Varying the 
normalization, a, and exponent, b, of the M-a relation such 
that a G [7,9] and b G [4,6] has very little impact on the envi¬ 
ronmental constraints. 

5.1.3. Constraints on SMBH binary population eccentricity 

It is not only the astrophysical environment of SMBH bi¬ 
naries that can induce a bend in the characteristic strain spec¬ 
trum. Binaries with non-zero eccentricity emit GWs at a spec¬ 
trum of harmonics of the orbital frequency. The cumulative 
effect over the entire population can lead to a depletion of 
the low frequency strain spectrum (Enoki et al. 2007; Sesana 
2013; Ravi et al. 2014; Huerta et al. 2015), and a turnover 
whose shape can be captured with the parametrized spectrum 
model employed in this paper. Hence, we can use our /tum 
posterior from the marginalization of /bend over all k, to de¬ 
duce constraints on the eccentricity of binaries at some refer¬ 
ence orbital separation. Our approach follows from the pre¬ 
vious astrophysics constraints, where we build populations 
and strain spectra which have varying eccentricities at a fixed 
semi-major axis of 0.01 pc, then construct a relationship be¬ 
tween this eccentricity and the spectral turnover. An impor¬ 
tant modeling assumption we make here is that binaries are 
(and always have been) driven entirely by GW emission. This 
factors into how we model d/^/dtr and how we evolve the 
binary eccentricity, where we assume binaries could have ec¬ 
centricities arbitrarily close to 1 in the past. 

We construct eccentric populations and corresponding 
strain spectra using the formalism of Huerta et al. (2015). 
The resulting relationship between the spectral turnover fre¬ 
quency and the eccentricity of all binaries at a semi-major 
axis of 0.01 pc is shown in the top panel of Fig. 11, along 



Figure 11. Same as Figure 10 except now we display the empirical map¬ 
ping (top) and posterior distribution (bottom) for the eccentricity of SMBH 
binaries when they had a semi-major axis of 0.01 pc. 

with the corresponding eccentricity posteriors from each am¬ 
plitude prior in the bottom panel. The high turnover frequency 
obtained with the MOP14 prior leads to an eccentricity pos¬ 
terior distribution that largely favors eo ^ 0.7 while the S13 
prior leads to an eccentricity posterior that is consistent with 
smaller eccentricities, more weakly favoring eo > 0.5. We 
emphasize that, whilst these eccentricities seem rather large, 
it is well established that binaries evolving in stellar environ¬ 
ments tend to increase their eccentricity (Quinlan 1996). It 
is therefore likely that most binaries can get to e ~ 0.5-0.7 
along their evolution (see tracks in Sesana 2010). The eccen¬ 
tricity growth rate is generally larger for smaller binary mass 
ratios, and for larger initial eccentricities. The latter is indeed 










NANOGrav Nine-year Isotropic GWB Limit 


17 


a crucial parameter; if, following galaxy mergers, the MBHB 
already has a significant (e > 0.5) eccentricity at the moment 
of formation, the subsequent evolution will almost certainly 
drive it to e > 0.9. Given that the MBHB eccentricity at for¬ 
mation is hard to determine (Aarseth 2003; Hemsendorf et al. 
2002; Berentzen et al. 2009; Amaro-Seoane et al. 2009), it is 
impossible to draw astrophysical conclusions from the con¬ 
strains above. Furthermore, in reality the history of a binary’s 
eccentricity will see phases of growth and circularization de¬ 
pending upon the interplay of environmental factors and GW 
emission (e.g., Sesana 2010; Kocsis & Sesana 2011). This 
should be considered in future analyses. 

5.2. Cosmological Model Limits 

SMBH binary mergers are not the only possible source of a 
GWB in the pulsar timing band. In this subsection, we use our 
power-law spectrum results, shown in Fig. 3, to constrain the 
fractional energy density of the Universe in relic GWs, along 
with stringent limits on the tension of cosmic strings. 

5.2.1. Relic GWs 

Quantum fluctuations of the gravitational field in the early 
universe, amplified by an inflationary phase, are expected 
to produce a stochastic relic GW background, see e.g. Gr- 
ishchuk (1976, 1977); Starobinsky (1980); Linde (1982). Ob¬ 
servations of this background would provide a unique in¬ 
sight into the highly curved spacetime of the early universe 
at less than 10“^^ s after the Big Bang and at energy scales 
of lO'® GeV, when quantum mechanics and general relativ¬ 
ity should reconcile, BICEP2/Keck et al. (2015); Ade et al. 
(2014). This background is expected to produce a character¬ 
istic signature in the polarization of the Cosmic Microwave 
Background (CMB) radiation, as well as CMB temperature 
anisotropies Grishchuk (2005). In the context of PTAs, the 
amplitude of the relic GW background is of astrophysical and 
cosmological interest due to the fact that it intrinsically de¬ 
pends on the equation of state of the early universe, w, and 
thus the Hubble constant in the inflationary stage //*, as well 
as the tensor-to-scalar ratio r , see e.g. Zhao (2011); Zhao et al. 
(2013). 

Specifically, the spectral index of the stochastic GWB, a, is 
related to the equation of state of the early universe, w, via 


where n, is the spectral index of the primordial power spec¬ 
trum, usually set to 0. In curi'ent hot big bang models, w = 1/3 
and rit = 0, thus a = -l. We therefore fix the spectral index to 
-1 and apply the Bayesian analysis methods described in Sec 
3.2 to the 9-year NANOGrav dataset. Using analysis meth¬ 
ods described in 4.2, we obtain a 95% upper limit of on the 
amplitude of the relic GW background of 

A/J/ = 8.1 X 10“'®, (37) 


Table 3 

Values for Qgy^{f)h^ reported at 1/yr. Cosmological parameters used for 
are h = 0.702, 7cmb = 0.276 K, = 0.725, = 0.275, and ^eq = 3454. 

Values of are in given in multiples of the Planck mass mpi = 1 / y/G. 


EoS, u) 

Spectral index, ol 

4 95% 



0 

-2 

1.2 X 10“** 

8.6 X 10“*2 

5.4 

1/3 

-1 

8.1 X 10“** 

4.2 X 10“*** 

1.6 X 10“2 

1 

-1/2 

2.0 X 10“** 

2.5 X 10““* 

8.2 X 10“^ 


where / is the frequency, pc = Stt/{3Hq) is the critical energy 
density required to close the Universe, Hq is the Hubble ex¬ 
pansion rate, and pg* is the total energy density in GWs (Allen 
& Romano 1999; Maggiore 2000). The NANOGrav limit is 
therefore 

f^gw(/)/r^ <4.2x 10“'°, (39) 

in a radiation-dominated universe with equation of state of w = 
1/3. This new limit is a factor of 2.9 better than the previous 
best limit of f2gw(/)/!^ = 1.2 x 10“® from LTM15. 

Although a radiation-dominated era is usually assumed to 
follow inflation in the hot big bang paradigm, there is cur¬ 
rently no evidence to show this held before Big Bang Nucle¬ 
osynthesis (BBN). In fact, the existence of a reheating stage or 
the existence of a cosmic phase transition both violate this as¬ 
sumption (Grishchuk 2001 ; Watanabe & Komatsu 2006; Zhao 
2011). For completeness, we now allow for other equations 
of state of the early universe before the BBN stage. The same 
analysis can be repeated assuming different equations of state 
for the early universe: a matter-dominated universe would 
have w = Q and therefore by Eq. (36), a = -2. If the universe 
were instead dominated by the kinetic energy of the inflaton, 
then w=l and a = -l /2. 

Einally, we place limits on the Hubble parameter, //* in the 
inflationary stage using methods developed in Zhao (2011). 
There, they introduced a way of translating the upper limit on 
a primordial GW background to a constraint on //». Using 
typical cosmological parameters h = 0.702, Tcmb = 0.276 K, 
Ha = 0.725, H,„ = 0.275, and Zeq = 3454, they obtain the fol¬ 
lowing relation: 

logioHgw(k*)= 1.25-:^2:^+21ogio ^ , (40) 

3w-l-t \fnp\ J 

where k* = 27r/* and is reported at a reference frequency of 
/» = 1 /yr and mpi = 1 / s/G is the Planck mass. Using Eq. 
(40), we can then limit for a fixed equation of state. Eor 
example, using the limit on = 4.2 x 10“'° for w = 1 /3, 
see Table 3, one can place a limit //* = 1.6 x 10“^mpi. Results 
for w = 0,1 are in Table 3. 

In LTM15, the limit on a primordial GW background with 
w = 1/3 and n, = 0 is Llg^(f)h^ = 1.2 x 10“°, resulting in = 
2.74 X 10“^mpi. One can see that the NANOGrav limit on 
//* is an improvement of 1.7 over the EPTA, and the most 
stringent limit to date. 


assuming a power spectrum for the characteristic strain with 
a = -1 at a reference frequency of / = yr“'. This then con¬ 
strains the GW energy density content per unit logarithmic 
frequency, divided by the critical energy density, pc, to close 
the Universe: 


f^gw(/) 


1 dpgw 
Pc dlog/ 


27r^ 

3H^ 


fhlif), 


(38) 


5.2.2. GW Background from Cosmic (Super)strings 

Cosmic strings are linear topological defects that can be 
produced in the early universe via phase transitions (Kib¬ 
ble 1976; Vilenkin 1981; Vilenkin & Shellard 2000). Cos¬ 
mic superstrings are fundamental strings stretched to cos¬ 
mological scales and arise in a wide range of string-theory- 
inspired cosmological scenarios (Sarangi &. Tye 2002; Jones 
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Figure 12. {left): Cosmic string constraints in terms of string tension Gfi and reconnection probability p using the results of recent cosmic string simulations 
described in Blanco-Pillado et al. (2014). (right): Cosmic string tension Gp vs loop size parameterized by acs using the model described in LTM15. The shaded 
area is ruled out by our GW upper limit in both panels. 


et al. 2003; Copeland et al. 2004). Cosmic (super)strings pro¬ 
duce a stochastic background of GWs as well as individual 
bursts (Damour & Vilenkin 2001; Damour & Vilenkin 2005; 
Siemens et al. 2006, 2007; Olmez et al. 2010). 

Our limits on the amplitude of the stochastic background 
can also be used to constrain the parameter space of cos¬ 
mic (super)strings. Recent simulations (Blanco-Pillado et al. 
2014) have shown that cosmic (super)string loop densities 
are dominated by loops that formed at scales comparable to 
the Hubble size at the time of formation, even though only 
about 10% of loops are formed with such large sizes. We use 
the loop distributions derived by Blanco-Pillado et al. (2014), 
specifically Eqs. (63), (65), and (67) of that reference with 
loop size acs = 0.05, together with the techniques described in 
Olmez et al. (2010) to compute the stochastic background pro¬ 
duced by cosmic string cusps. The cosmological parameters 
we used are taken from the Planck 2015 data (Ade et al. 2015). 
In this case the relevant parameters are the string tension Gp 
and the reconnection probability p. We explore this parame¬ 
ter space and exclude regions where the cosmic (super)string 
network would have resulted in a stochastic background am¬ 
plitude larger than that ruled out by our measurements. The 
left panel of Figure 12 shows the results of our analysis. On 
the y-axis we show the reconnection probability and on the x- 
axis the string tension. The gray shaded area shows the region 
of cosmic string parameter space that is ruled out. Note that 
for p = 1 our data only allow for cosmic (super)strings with 
tensions Gp < 1.3 x 10“'°. 

Recently LTM15 presented a comprehensive and fully gen¬ 
eral overview of cosmic string limits from the EPTA, and 
found a conservative limit on the string tension to be G/i < 
1.3 X 10“^. The limit is conservative in the sense that it is 
found by considering a wide range of loop sizes and taking 
the upper limit to be the largest possible value of G/i con¬ 
sistent with the data. The limit was identical to that set by 
the Planck Collaboration, combining Planck and high-/ Cos¬ 
mic Microwave Background data with Atacama Cosmology 
Telescope (ACT) and the South Pole Telescope (SPT) , cf. 
Planck Collaboration et al. (2014). While the calculation in 
LTM15 was not carried out explicitly for the Blanco-Pillado 
et al. (2014) simulations we can use their published limit 
on f2gw(/)/!^ = 1.2 X 10 ° for cosmic strings to place a limit 
of Gfi < 8.6 X 10 '°. Our limit for this model is therefore 


roughly a factor of 6.6 times more constraining than the in¬ 
ferred previous limit. Using the same analysis developed by 
the EPTA, (Sanidas et al. 2012, 2013), we compute the up¬ 
per limit on the string tension Gp, as a function of loop size 
Ucs as shown in the right panel of Figure 12. Our conserva¬ 
tive limit on cosmic string tension using this range of cosmic 
string models is Gp < 3.3 x 10“*', a factor of 4 better than both 
the combined Planck, ACT, SPT limit and the EPTA limit. 

6. SUMMARY AND CONCLUSIONS 

This paper reports on the search for an isotropic stochas¬ 
tic GW background in NANOGrav’s 9-year dataset. We do 
not find positive statistical evidence for the presence of such 
a signal. Following up on a series of earlier results by the 
three PTAs, we report new upper limits on the amplitude of 
backgrounds described by power-law spectra; 

• For an astrophysical background of SMBH binaries 
(corresponding to a timing-residual spectral density 
with exponent 7 = 13/3), we find a 95% confidence 
limit Agw < 1.5 X 10“*^, five times more constrain¬ 
ing than the analogous limit for NANOGrav’s 5-year 
dataset (DFG13). Under the assumption of purely 
GW-driven evolution, leading to an unbroken 7=13/3 
power law, we compute the probability that our con¬ 
straint is consistent with the MOP14 and S13/RWS14 
theoretical predictions for Agw as 0 . 8 % and 20 %, re¬ 
spectively, essentially ruling out the MOP14 model and 
placing the S13/RWS14 model in tension with our data. 
[Sec. 4.2.1.] 

• We verify the consistency of our limit with previously 
reported scaling relations between SMBH mass and 
galactic bulge mass, adopting fiducial estimates for 
galaxy merger rates and the stellar mass function. Un¬ 
der the assumption of circular GW-driven binaries, our 
limit is slightly inconsistent with the Kormendy & Ho 
(2013) relation, and consistent within the error margin 
for the McConnell & Ma (2013) relation. [Sec. 5.1.1.] 

• We also perform an optimal-statistic (cross-correlation) 
analysis, and find limits that are 5.4 and 1.5 times more 
constraining than the analogous DEG 13 and LTM15 re¬ 
sults. The cross-correlation SNR is 0.8, indicating that 
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there is little evidence for inter-pulsar residual correla¬ 
tions induced by GWs. [Sec. 4.1.] 

• We extend the power-law background search to generic 
spectral indices, and place the most stringent limits 
so far on the energy density of relic GWs, flgv/if = 
l/yr)/;^ < 4.2 x 10“'°, for a w = 1/3 early-Universe 
equation of state. From this we obtain limits on the 
Hubble parameter during inflation, //* = 1.6 X 10 ^ nipi. 
[Secs. 4.2.1, 5.2.1] 

• We place the most stringent limits to date on a GW 
background generated by a network of cosmic strings, 
^gwif = l/yr)/i^ < 2.2 x 10“'°, which translates into 
a conservative upper limit on cosmic string tension 
Gfi < 3.3 X 10“*, using the model presented in LTM15. 
This is a factor of 4 better than both the combined 
Planck, ACT, SPT limit and the EPTA limit. Using the 
recent models of Blanco-Pillado et al. (2014) we find 
G/i < 1.3 X 10“'°, a factor of 6.6 times more constrain¬ 
ing than an identical analysis using the EPTA limit. 
[Secs. 4.2.1, and 5.2.2.] 

We further probe the interface between PTA observations 
and SMBH-binary population estimates by analyzing the 9- 
year dataset in terms of a GW background described by an 
inflected power law [Eq. (24)]. We derive joint posteriors for 
the spectral parameters (the amplitude Agw, the inflection fre¬ 
quency /bend, and the shape exponent 7) assuming Ag* priors 
from MOP14 and S13/RWS14 (see Eig. 5). Eor both priors 
(but more so for MOP14), the inflected power-law model is 
prefetTed to an unbroken power law. The /bend posterior can 
be used to infer astrophysical information about the effects 
that may reduce GW power at lower frequencies, such as the 
initial eccentricity of SMBH binaries, and the environmental 
influences of stars and gas in galactic nuclei. To summarize: 

• We find that the data prefer an inflected spectrum over 
a pure power law, with Bayes factors ~ 22 and 2.2 
for the MOP14 and S13 amplitude priors, respectively. 
The same analysis, run on the 5-year DEG13 dataset, 
provides little to no information about the shape of the 
GWB spectrum. [Sec. 4.2.1.] 

• Under several simplifying assumptions, we map the 
posterior distribution of /bend into posterior distribu¬ 
tions for the nuclear stellar mass density p (which mod¬ 
ulates the strength of binary frequency evolution by 
stellar scattering, cf Sec. 5.1.2), the SMBH mass ac¬ 
cretion rate from circumbinary disks M\ (which can 
be linked to binary frequency evolution by interactions 
with gas, cf. Sec. 5.1.2), and the initial value of SMBH- 
binary orbital eccentricity ea (which distributes GW 
power to higher frequency harmonics, cf. Sec. 5.1.3). 

Eor the last decade (and longer), the three PTA consortia 
have been engaged in an accelerating race toward higher GW 
sensitivities—the analysis presented in this paper represents 
the latest stage of the race, but not the last. While our investi¬ 
gation cannot claim the ultimate prize, a positive detection, it 
is the first to use information about the amplitude and shape 
of GW background to make concrete (if limited) astrophysical 
statements about the dynamics and environments of SMBH 
binaries. The era of GW astrophysics is truly upon us. 
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